%reads the state of model ./outfiles/-1.txt in terms of the marker
%properties: temperature, stress and strain rate

%input: positions of new markers, mx,mz and position of the new grid
%output: marker temperatures mTemp, stresses, mstressxx, mstressxz and
%strain rates, mstrainxx,mstrainxz and a surface with the position of the
%LAB,xLAB,zLAB

%LAB is defined as the 1200 degC isotherm and xLAB is from the grid of
%-1.txt

%omgrid makenewmodel, hvis problemer

function [mTemp,mstressxx,mstressxz,mstrainxx,mstrainxz,xLAB,zLAB]=importstate(mx,mz,xmax,zmax);
N=length(mx);

%loading from -1.txt
count=-1;
dxold=readdata(count,8);
dzold=readdata(count,9);
mxold=readdata(count,10);
mzold=readdata(count,11);
mTemp=readdata(count,15);

Sold=length(dxold)+1;
xold=zeros(1,Sold);
for s=1:Sold-1 xold(s+1)=xold(s)+dxold(s); end

%if max(xold)~=xmax ....

[grid,xgrid,zgrid]=bilingrid(mTemp,mxold,mzold,dxold,dzold,1,1); [Zgrid,Xgrid] = meshgrid(zgrid,xgrid);
mTemp=interp2(Xgrid',Zgrid',grid',mx,mz); mTemp(find(isnan(mTemp)))=max(mTemp);
cont=contourc(xgrid,zgrid,grid',[1200 1200]); %here the LAB temp of 1200 degC is defined
pos=1; xcont=[]; ycont=[];
while pos<length(cont)
    step=cont(2,pos);
    xcont=[xcont cont(1,pos+1:pos+step)]; ycont=[ycont cont(2,pos+1:pos+step)];
    pos=pos+step+1;
end
[xcont sortedI]=sort(xcont); ycont=ycont(sortedI); 
xcont=[0 xcont xmax]; ycont=[ycont(1) ycont ycont(end)];
[xcont sortedI]=unique(xcont); ycont=ycont(sortedI); %removing duplicant x vals
zLAB=interp1(xcont,ycont,xold); xLAB=xold;
figure; contour(Xgrid,Zgrid,grid,100);  set(gca,'yDir','rev'); shading flat; hold on; plot(xold,zLAB,'k','linewidth',3); 

mstressxx=readdata(count,16); [grid,xgrid,zgrid]=bilingrid(mstressxx,mxold,mzold,dxold,dzold,1,1); [Z,X] = meshgrid(zgrid,xgrid);
mstressxx=interp2(X',Z',grid',mx,mz); mstressxx(find(isnan(mstressxx)))=1;
mstressxz=readdata(count,17); [grid,xgrid,zgrid]=bilingrid(mstressxz,mxold,mzold,dxold,dzold,1,1); [Z,X] = meshgrid(zgrid,xgrid);
mstressxz=interp2(X',Z',grid',mx,mz); mstressxz(find(isnan(mstressxz)))=0;

mstrainxx=readdata(count,18); [grid,xgrid,zgrid]=bilingrid(mstrainxx,mxold,mzold,dxold,dzold,1,1); [Z,X] = meshgrid(zgrid,xgrid);
mstrainxx=interp2(X',Z',grid',mx,mz); mstrainxx(find(isnan(mstrainxx)))=0;
mstrainxz=readdata(count,19); [grid,xgrid,zgrid]=bilingrid(mstrainxz,mxold,mzold,dxold,dzold,1,1); [Z,X] = meshgrid(zgrid,xgrid);
mstrainxz=interp2(X',Z',grid',mx,mz); mstrainxz(find(isnan(mstrainxz)))=0;



% S=201; T=201; dim=[S T];
% xmax=max(xout); zmax=max(mz); N=length(mx);
% dx=xmax/(S-1); dz=zmax/(T-1); x=0:dx:xmax; z=0:dz:zmax;
% transitionwidth=5e4;
% rhouc=2800; rholc=3000; rhom=3300; alpha=0.5; %the fraction of upper to lower crust
% thinliththickness=100e3; thickliththickness=-50e3; meanthickness=0.5*(thinliththickness+thickliththickness);
% %liththickness=thinliththickness+thickliththickness*exp(-(x-xmax/2).^2./(2*transitionwidth^2));
% %x1=9e5; x2=11e5; liththickness=150e3-0.5*0*50e3*(erf(1/transitionwidth*(x-x1))-erf(1/transitionwidth*(x-x2))+2);
% liththickness=150e3+1e3*sin(2*pi/3e5*x);
% %liththickness=meanthickness+(meanthickness-thinliththickness)*erf(1/transitionwidth*(x-xmax/2));
% thincrustthickness=30e3; thickcrustthickness=40e3; meanthickness=0.5*(thincrustthickness+thickcrustthickness);
% crustthickness=meanthickness+(meanthickness-thincrustthickness)*erf(1/transitionwidth*(x-xmax/2));
% %crustthickness=thincrustthickness+thickcrustthickness*exp(-(x-xmax/2).^2./(2*transitionwidth^2));
% 
% Tsurface=0;
% 
% C1=(alpha*rhouc+(1-alpha)*rholc)*thincrustthickness;
% C2=(alpha*rhouc+(1-alpha)*rholc)*crustthickness;
% uplift=(C1-C2+(crustthickness-thincrustthickness)*rhom)/rhom;
% weaklayerthickness=40000-uplift;
% 
% %midcrustdepth=0.5*alpha*crustthickness.*(1-erf(1/transitionwidth*(x-xmax/2)))+weaklayerthickness;
% midcrustdepth=0.5*alpha*thincrustthickness*(erf(1/transitionwidth*(x-900e3))+1)-0.5*alpha*thincrustthickness*(erf(1/transitionwidth*(x-1100e3))+1)+weaklayerthickness;
% mohodepth=crustthickness+weaklayerthickness;
% lithdepth=liththickness+weaklayerthickness;
% 
% 
% % for s=2:floor(S/2)
% %     temp(s,:)=temp(1,:);
% % end
% % for s=floor(S/2)+1:S-1
% %     temp(s,:)=temp(S,:);
% % end
% 
% count=-1;
% 
% dxold=readdata(count,8);
% dzold=readdata(count,9);
% mxold=readdata(count,10);
% mzold=readdata(count,11);
% mTemp=readdata(count,15);
% 
% [grid,xgrid,zgrid]=bilingrid(mTemp,mxold,mzold,dxold,dzold,1,1); [Z,X] = meshgrid(zgrid,xgrid);
% mTemp=interp2(X',Z',grid',mx,mz); mTemp(find(isnan(mTemp)))=max(mTemp);
% cont=contourc(xgrid,zgrid,grid',[1200 1200]);
% pos=1; xcont=[]; ycont=[];
% while pos<length(cont)
%     step=cont(2,pos);
%     xcont=[xcont cont(1,pos+1:pos+step)]; ycont=[ycont cont(2,pos+1:pos+step)];
%     pos=pos+step+1;
% end
% [xcont sortedI]=sort(xcont); ycont=ycont(sortedI); 
% xcont=[0 xcont xmax]; ycont=[ycont(1) ycont ycont(end)];
% [xcont sortedI]=unique(xcont); ycont=ycont(sortedI); %removing duplicant x vals
% lithdepth=interp1(xcont,ycont,x);
% figure(1); contour(X,Z,grid,100);  set(gca,'yDir','rev'); shading flat; hold on; plot(x,lithdepth,'k','linewidth',3); plot(x,mohodepth,'k','linewidth',3); plot(x,midcrustdepth,'k','linewidth',3); plot(x,weaklayerthickness,'k','linewidth',3); colorbar;
% 
% mstressxx=readdata(count,16); [grid,xgrid,zgrid]=bilingrid(mstressxx,mxold,mzold,dxold,dzold,1,1); [Z,X] = meshgrid(zgrid,xgrid);
% mstressxx=interp2(X',Z',grid',mx,mz); mstressxx(find(isnan(mstressxx)))=1;
% mstressxz=readdata(count,17); [grid,xgrid,zgrid]=bilingrid(mstressxz,mxold,mzold,dxold,dzold,1,1); [Z,X] = meshgrid(zgrid,xgrid);
% mstressxz=interp2(X',Z',grid',mx,mz); mstressxz(find(isnan(mstressxz)))=0;
% 
% mstrainxx=readdata(count,18); [grid,xgrid,zgrid]=bilingrid(mstrainxx,mxold,mzold,dxold,dzold,1,1); [Z,X] = meshgrid(zgrid,xgrid);
% mstrainxx=interp2(X',Z',grid',mx,mz); mstrainxx(find(isnan(mstrainxx)))=0;
% mstrainxz=readdata(count,19); [grid,xgrid,zgrid]=bilingrid(mstrainxz,mxold,mzold,dxold,dzold,1,1); [Z,X] = meshgrid(zgrid,xgrid);
% mstrainxz=interp2(X',Z',grid',mx,mz); mstrainxz(find(isnan(mstrainxz)))=0;
% 
% %index=find((mx-xmax/2).^2+(mz-500e3).^2<(50e3)^2); mTemp(index)=1600;
% 
% mtype=0*mTemp+3; %3 is the mantle
% Iweaklayer=find(interp1(x,weaklayerthickness,mx)>mz); mtype(Iweaklayer)=0; %0 is the weak layer
% Iucrust=find(interp1(x,weaklayerthickness,mx)<=mz & interp1(x,midcrustdepth,mx)>mz); mtype(Iucrust)=1; %1 is the upper crust
% Ilcrust=find(interp1(x,midcrustdepth,mx)<=mz & interp1(x,mohodepth,mx)>mz); mtype(Ilcrust)=2; %2 is the lower crust
% %Iload=find(abs(mx-xmax/2)<50e3 & mz>interp1(x,weaklayerthickness,mx)-10e3 & mz<=interp1(x,weaklayerthickness,mx)); mtype(Iload)=5;
% loadthickness=5e3*exp(-(x-xmax/2).^2./(2*100e3^2));
% %Iload=find(abs(mz>interp1(x,weaklayerthickness-loadthickness,mx)) & mz<=interp1(x,weaklayerthickness,mx)); mtype(Iload)=5;
% Iast=find(interp1(x,lithdepth,mx)<mz); mtype(Iast)=4; %4 is the asthenosphere
% 
% Irightlith=find(mtype==3 & mx<xmax/2); mtype(Irightlith)=4;
% 
% %Iload=find(abs(mx-xmax/2)<50e3 & mz>interp1(x,weaklayerthickness,mx)-20e3 & mz<=interp1(x,mohodepth,mx)); mtype(Iload)=1;
% %Iload=find(abs(mx-xmax/2)<50e3 & mz>interp1(x,weaklayerthickness,mx)+100e3 & mz<=interp1(x,mohodepth+100e3,mx)); mtype(Iload)=1;
% %Iload=find(abs(mx-xmax/2)<25e3 & abs(mz-zmax/2)<25e3); mtype(Iload)=1;
% 
% 
% surface=interp1(x,weaklayerthickness,xout);
% midcrust=interp1(x,midcrustdepth,xout);
% moho=interp1(x,mohodepth,xout);
% topast=interp1(x,lithdepth,xout);

end